# Summary: This script replicates the article's main analysis

#####################################################
#----------Replicate Main Paper Results ------------#
#####################################################

options(scipen=999)
rm(list = ls())

# Load libraries
library(rddensity)
library(rdrobust)
library(rdpower)
library(sandwich)
library(stargazer)
library(readxl)
library(lmtest)

# Set WD
# setwd('~/Dataverse/')

# Source functions
source('./Code/fun_analysis.R')

# Load Data
load(file = "./Data/d.RData")

# Keep only observations from municipalities where both parties had male and female candidates
d1 <- subset(d, nonZeroDF == 2)
length(unique(d1$codeyear))

# Get codeyear for where a centralized and decentralized party won
centParties <- d1$codeyear[(d1$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d1$eleito == 1]
decentParties <- d1$codeyear[!(d1$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d1$eleito == 1]

#################################################################
####### Only Centralized Parties are able to deliver votes ######
#################################################################

# FIGURE 1
pdf('./MainPaperResults/fig1.pdf')
# All Parties
genPlot(mydata = subset(d1), 
	DV = 'avg_DF', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1)$codeyear)
# Centralized
genPlot(mydata = subset(d1, codeyear %in% centParties), 
	DV = 'avg_DF', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1, codeyear %in% centParties)$codeyear)
# Descentralized
genPlot(mydata = subset(d1, codeyear %in% decentParties), 
	DV = 'avg_DF', 
	dvName = 'Federal Deputy', 
	vs = T, 
	clusterVar = subset(d1, codeyear %in% decentParties)$codeyear)
dev.off()

# Main Models: h = optimal 
out.cent.opt <- rdrobust(y = subset(d1, codeyear %in% centParties)$avg_DF
		, x = subset(d1, codeyear %in% centParties)$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties)$codeyear
		, all = TRUE)
out.decent.opt <- rdrobust(y = subset(d1, codeyear %in% decentParties)$avg_DF
		, x = subset(d1, codeyear %in% decentParties)$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties)$codeyear
		, all = TRUE)
# Test Difference:
zOut <- zDif(beta1 = out.cent.opt$coef[3,], beta2 = out.decent.opt$coef[3,]
	, se1 = out.cent.opt$se[3,], se2 = out.decent.opt$se[3,])
zOut
(1 - pnorm(zOut$zscore)) * 2

#### Results by Period
d1$period <- 'weak'
d1$period[d1$ANO_ELEICAO > 2007] <- 'strong'
d1$period[d1$ANO_ELEICAO > 2015] <- 'moderate'

# Main Models: h = optimal 
out.cent.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_DF
		, x = subset(d1, codeyear %in% centParties & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)
out.cent.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_DF
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)
out.cent.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_DF
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)
out.decent.opt.weak <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'weak')$avg_DF
		, x = subset(d1, codeyear %in% decentParties & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'weak')$codeyear
		, all = TRUE)
out.decent.opt.strong <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'strong')$avg_DF
		, x = subset(d1, codeyear %in% decentParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'strong')$codeyear
		, all = TRUE)
out.decent.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% decentParties & period == 'moderate')$avg_DF
		, x = subset(d1, codeyear %in% decentParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% decentParties & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_cen <- c(out.cent.opt.weak$coef[3,], out.cent.opt.strong$coef[3,]
	, out.cent.opt.moderate$coef[3,])
se_cen <- c(out.cent.opt.weak$se[3,], out.cent.opt.strong$se[3,]
	, out.cent.opt.moderate$se[3,])
lb_cen <- c(out.cent.opt.weak$ci[3,1], out.cent.opt.strong$ci[3,1]
	, out.cent.opt.moderate$ci[3,1])
ub_cen <- c(out.cent.opt.weak$ci[3,2], out.cent.opt.strong$ci[3,2]
	, out.cent.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_decen <- c(out.decent.opt.weak$coef[3,], out.decent.opt.strong$coef[3,]
	, out.decent.opt.moderate$coef[3,])
se_decen <- c(out.decent.opt.weak$se[3,], out.decent.opt.strong$se[3,]
	, out.decent.opt.moderate$se[3,])
lb_decen <- c(out.decent.opt.weak$ci[3,1], out.decent.opt.strong$ci[3,1]
	, out.decent.opt.moderate$ci[3,1])
ub_decen <- c(out.decent.opt.weak$ci[3,2], out.decent.opt.strong$ci[3,2]
	, out.decent.opt.moderate$ci[3,2])

# Calculate Differences 
zOutDecent <- zDif(beta1 = out.decent.opt.strong$coef[3,], beta2 = out.decent.opt.weak$coef[3,]
	, se1 = out.decent.opt.strong$se[3,], se2 = out.decent.opt.weak$se[3,])
zOutCent <- zDif(beta1 = out.cent.opt.strong$coef[3,], beta2 = out.cent.opt.weak$coef[3,]
	, se1 = out.cent.opt.strong$se[3,], se2 = out.cent.opt.weak$se[3,])


# Start the plot:
# FIGURE 2
pdf('./MainPaperResults/fig2.pdf', width = 10)
par(mar = c(5, 5, 2, 2))
plot(x = NULL, y = NULL, axes = FALSE, 
	cex = 2, xlim = c(0.75, 2.25), ylim = c(-2, 2),
	ylab = "RD Effect of Co-Partisan Mayor", col = c('black', 'grey46'), 
	pch = c(20, 17),
	xlab = "Party Type", 
	cex.lab = 2)
sapply(1:2, function(i) lines(x = c(c(0.9, 1.1)[i], c(0.9, 1.1)[i]), 
	y = c(lb_decen[i], ub_decen[i]), col = c('black', 'gray46')[i]))
sapply(1:2, function(i) lines(x = c(c(1.9, 2.1)[i], c(1.9, 2.1)[i]), 
	y = c(lb_cen[i], ub_cen[i]), col = c('black', 'gray46')[i]))
axis(1, at = seq(1, 2, 1), label = c('Decentralized', 'Centralized'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_decen[1], y1 = betas_decen[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_decen[1], y1 = betas_decen[2], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_decen[2], y1 = betas_decen[2], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_cen[1], y1 = betas_cen[1], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_cen[1], y1 = betas_cen[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_cen[2], y1 = betas_cen[2], lty = 2)

# Add points for centralized
points(x = c(1.9, 2.1), y = betas_cen[1:2], pch = c(20, 17), col = c('black', 'grey46'), cex = 2)
# Add points for decentralized
points(x = c(0.9, 1.1), y = betas_decen[1:2], pch = c(20, 17), col = c('black', 'grey46'), cex = 2)


# Add information about difference
text(round(zOutDecent$difference, 3), y = betas_decen[1] + 0.75, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutDecent$zscore > 0, (1 - pnorm(zOutDecent$zscore)) * 2, pnorm(zOutDecent$zscore) * 2), 3), ')')
, y = betas_decen[1] + 0.5, x = 1, cex = 1.5)

text(round(zOutCent$difference, 3), y = betas_cen[1] + 1.15, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutCent$zscore > 0, (1 - pnorm(zOutCent$zscore)) * 2, pnorm(zOutCent$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutCent$zscore > 0, (1 - pnorm(zOutCent$zscore)) * 2, pnorm(zOutCent$zscore) * 2), 3)), ')')
, y = betas_cen[1] + 0.9, x = 2, cex = 1.5)

# Add legend
legend('bottomleft', 
	title = 'Electoral Rule Regime:',
	legend = c('Strong (2010-2014)', 'Weak (1998-2006)'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()


##################################################################################
####### H1: male candidates > female candidates By Period, only centralized ######
##################################################################################

d1$period <- 'weak'
d1$period[d1$ANO_ELEICAO > 2007] <- 'strong'
d1$period[d1$ANO_ELEICAO > 2015] <- 'moderate'

# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_male_DF
		, x = subset(d1, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'weak')$avg_fem_DF
		, x = subset(d1, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'weak')$codeyear
		, all = TRUE)

# Strong
out.cent.men.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_male_DF
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'strong')$avg_fem_DF
		, x = subset(d1, codeyear %in% centParties & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'strong')$codeyear
		, all = TRUE)

# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_male_DF
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% centParties & period == 'moderate')$avg_fem_DF
		, x = subset(d1, codeyear %in% centParties & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d1, codeyear %in% centParties & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
	, out.cent.men.opt.moderate$coef[3,])
se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
	, out.cent.men.opt.moderate$se[3,])
lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
	, out.cent.men.opt.moderate$ci[3,1])
ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
	, out.cent.men.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
	, out.cent.women.opt.moderate$coef[3,])
se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
	, out.cent.women.opt.moderate$se[3,])
lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
	, out.cent.women.opt.moderate$ci[3,1])
ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
	, out.cent.women.opt.moderate$ci[3,2])

# Calculate Differences (Paper's TABLE 1)
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
	, se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
	, se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
	, se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])

sink("./MainPaperResults/table1.txt")
print("TABLE 1")
strong_weak <- zDif(beta1 = zOutStrong$difference, se1 = zOutStrong$SE,
	beta2 = zOutWeak$difference, se2 = zOutWeak$SE)
strong_weak
(1 - pnorm(strong_weak$zscore))

strong_moderate <- zDif(beta1 = zOutStrong$difference, se1 = zOutStrong$SE,
	beta2 = zOutModerate$difference, se2 = zOutModerate$SE)
strong_moderate
(1-pnorm(strong_moderate$zscore))

moderate_weak <- zDif(beta2 = zOutWeak$difference, se2 = zOutWeak$SE,
	beta1 = zOutModerate$difference, se1 = zOutModerate$SE)
moderate_weak
(1 - pnorm(moderate_weak$zscore))
sink()

# Start the plot:
# FIGURE 3
pdf('./MainPaperResults/fig3.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
	pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-3, 3),
	ylab = "RD Effect of Co-Partisan Mayor",
	xlab = "Electoral Rule Regime", 
	cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
	y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
	y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = betas_men[1] + 0.5, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
, y = betas_men[1] + 0.2, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = betas_men[2] + 0.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
, y = betas_men[2] + 0.45, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = betas_men[3] + 0.5, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
, y = betas_men[3] + 0.2, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -1.5, 
	legend = c('Women', 'Men'), 
	col = c('gray46', 'black'), 
	pch = c(17, 20), 
	lty = c(1, 1),
	bty = 'n',
	cex = 1.5)
dev.off()


#####################################################################################
####### H2: male brokers > female brokers on Male VS, only centralized parties#######
#####################################################################################

# Select mixed-gender races
d$female <- as.numeric(d$CODIGO_SEXO == 4)
mixed <- aggregate(d$female, list(d$codeyear), sum)
names(mixed) <- c('codeyear', 'mixed')
mixed$mixed <- ifelse(mixed$mixed == 1, 1, 0)
d <- merge(d, mixed, by = 'codeyear')
# Subset to Mixed Races and parties ran both male and female candidates
d2 <- subset(d, mixed == 1 & nonZeroDF == 2)
# Create running variable
d2$run_var_fem <- ifelse(d2$position == 1 & d2$CODIGO_SEXO == 4
                         , d2$diff_cand * 1, d2$diff_cand * -1)
d2$run_var_fem[d2$mayorFem == 0] <- NA
d2$run_var_male <- ifelse(d2$position == 1 & d2$CODIGO_SEXO == 2
                          , d2$diff_cand * 1, d2$diff_cand * -1)
d2$run_var_male[d2$mayorFem == 1] <- NA
table(is.na(d2$run_var_male))
table(is.na(d2$run_var_fem))

#### Results by Period
d2$period <- 'weak'
d2$period[d2$ANO_ELEICAO > 2007] <- 'strong'
d2$period[d2$ANO_ELEICAO > 2015] <- 'moderate'

##############################
######### mayor = men ########
##############################
# Weak
out.h2.men.men.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'weak')$avg_male_DF
                                   , x = subset(d2, codeyear %in% centParties  & period == 'weak')$run_var_male
                                   , p = 1
                                   , c = 0	
                                   , cluster = subset(d2, codeyear %in% centParties & period == 'weak')$codeyear
                                   , all = TRUE)
out.h2.men.women.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'weak')$avg_fem_DF
                                     , x = subset(d2, codeyear %in% centParties & period == 'weak')$run_var_male
                                     , p = 1
                                     , c = 0	
                                     , cluster = subset(d2, codeyear %in% centParties & period == 'weak')$codeyear
                                     , all = TRUE)

# Difference 
outZ <- zDif(beta1 = out.h2.men.men.non.opt$coef[3,], beta2 = out.h2.men.women.non.opt$coef[3,]
             , se1 = out.h2.men.men.non.opt$se[3,], se2 = out.h2.men.women.non.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2
summary(out.h2.men.men.non.opt)
summary(out.h2.men.women.non.opt)

# Strong
out.h2.men.men.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'strong')$avg_male_DF
                                    , x = subset(d2, codeyear %in% centParties & period == 'strong')$run_var_male
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, codeyear %in% centParties  & period == 'strong')$codeyear
                                    , all = TRUE)
out.h2.men.women.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'strong')$avg_fem_DF
                                      , x = subset(d2, codeyear %in% centParties & period == 'strong')$run_var_male
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d2, codeyear %in% centParties & period == 'strong')$codeyear
                                      , all = TRUE)
summary(out.h2.men.men.with.opt)
summary(out.h2.men.women.with.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.men.men.with.opt$coef[3,], beta2 = out.h2.men.women.with.opt$coef[3,]
             , se1 = out.h2.men.men.with.opt$se[3,], se2 = out.h2.men.women.with.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2


# Moderate
out.h2.men.men.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'moderate')$avg_male_DF
                                       , x = subset(d2, codeyear %in% centParties & period == 'moderate')$run_var_male
                                       , p = 1
                                       , c = 0	
                                       , cluster = subset(d2, codeyear %in% centParties & period == 'moderate')$codeyear
                                       , all = TRUE)
out.h2.men.women.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'moderate')$avg_fem_DF
                                         , x = subset(d2, codeyear %in% centParties & period == 'moderate')$run_var_male
                                         , p = 1
                                         , c = 0	
                                         , cluster = subset(d2, codeyear %in% centParties & period == 'moderate')$codeyear
                                         , all = TRUE)
summary(out.h2.men.men.withLeg.opt)
summary(out.h2.men.women.withLeg.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.men.men.withLeg.opt$coef[3,], beta2 = out.h2.men.women.withLeg.opt$coef[3,]
             , se1 = out.h2.men.men.withLeg.opt$se[3,], se2 = out.h2.men.women.withLeg.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2

# Get point estimates, se, and ci for men
betas_men <- c(out.h2.men.men.non.opt$coef[3,], out.h2.men.men.with.opt$coef[3,]
               , out.h2.men.men.withLeg.opt$coef[3,])
se_men <- c(out.h2.men.men.non.opt$se[3,], out.h2.men.men.with.opt$se[3,]
            , out.h2.men.men.withLeg.opt$se[3,])
lb_men <- c(out.h2.men.men.non.opt$ci[3,1], out.h2.men.men.with.opt$ci[3,1]
            , out.h2.men.men.withLeg.opt$ci[3,1])
ub_men <- c(out.h2.men.men.non.opt$ci[3,2], out.h2.men.men.with.opt$ci[3,2]
            , out.h2.men.men.withLeg.opt$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.h2.men.women.non.opt$coef[3,], out.h2.men.women.with.opt$coef[3,]
                 , out.h2.men.women.withLeg.opt$coef[3,])
se_women <- c(out.h2.men.women.non.opt$se[3,], out.h2.men.women.with.opt$se[3,]
              , out.h2.men.women.withLeg.opt$se[3,])
lb_women <- c(out.h2.men.women.non.opt$ci[3,1], out.h2.men.women.with.opt$ci[3,1]
              , out.h2.men.women.withLeg.opt$ci[3,1])
ub_women <- c(out.h2.men.women.non.opt$ci[3,2], out.h2.men.women.with.opt$ci[3,2]
              , out.h2.men.women.withLeg.opt$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.h2.men.men.non.opt$coef[3,], beta2 = out.h2.men.women.non.opt$coef[3,]
                 , se1 = out.h2.men.men.non.opt$se[3,], se2 = out.h2.men.women.non.opt$se[3,])
zOutStrong <- zDif(beta1 = out.h2.men.men.with.opt$coef[3,], beta2 = out.h2.men.women.with.opt$coef[3,]
                   , se1 = out.h2.men.men.with.opt$se[3,], se2 = out.h2.men.women.with.opt$se[3,])
zOutModerate <- zDif(beta1 = out.h2.men.men.withLeg.opt$coef[3,], beta2 = out.h2.men.women.withLeg.opt$coef[3,]
                     , se1 = out.h2.men.men.withLeg.opt$se[3,], se2 = out.h2.men.women.withLeg.opt$se[3,])

# Start the plot:
# FIGURE 4.a
pdf('./MainPaperResults/fig4_a.pdf', width = 8)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
     pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-10, 10),
     ylab = "RD Effect of Co-Partisan Mayor",
     xlab = "Electoral Rule Regime", 
     cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
                              y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
                              y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = 7.75, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
     , y = 6.5, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = 7.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
     , y = 6.5, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = 7.75, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
     , y = 6.5, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -4.5, 
       legend = c('Women', 'Men'), 
       col = c('gray46', 'black'), 
       pch = c(17, 20), 
       lty = c(1, 1),
       bty = 'n',
       cex = 1.5)
dev.off()

################################
######### mayor = women ########
################################
# Weak
out.h2.women.men.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'weak')$avg_male_DF
                                     , x = subset(d2, codeyear %in% centParties & period == 'weak')$run_var_fem
                                     , p = 1
                                     , c = 0	
                                     , cluster = subset(d2, codeyear %in% centParties & period == 'weak')$codeyear
                                     , all = TRUE)
out.h2.women.women.non.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'weak')$avg_fem_DF
                                       , x = subset(d2, codeyear %in% centParties & period == 'weak')$run_var_fem
                                       , p = 1
                                       , c = 0	
                                       , cluster = subset(d2, codeyear %in% centParties & period == 'weak')$codeyear
                                       , all = TRUE)
summary(out.h2.women.men.non.opt)
summary(out.h2.women.women.non.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.non.opt$coef[3,], beta2 = out.h2.women.women.non.opt$coef[3,]
             , se1 = out.h2.women.men.non.opt$se[3,], se2 = out.h2.women.women.non.opt$se[3,])
outZ
(pnorm(outZ$zscore)) * 2

# Strong
out.h2.women.men.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'strong')$avg_male_DF
                                      , x = subset(d2, codeyear %in% centParties & period == 'strong')$run_var_fem
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d2, codeyear %in% centParties & period == 'strong')$codeyear
                                      , all = TRUE)
out.h2.women.women.with.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'strong')$avg_fem_DF
                                        , x = subset(d2, codeyear %in% centParties & period == 'strong')$run_var_fem
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(d2, codeyear %in% centParties & period == 'strong')$codeyear
                                        , all = TRUE)
summary(out.h2.women.men.with.opt)
summary(out.h2.women.women.with.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.with.opt$coef[3,], beta2 = out.h2.women.women.with.opt$coef[3,]
             , se1 = out.h2.women.men.with.opt$se[3,], se2 = out.h2.women.women.with.opt$se[3,])
outZ
(pnorm(outZ$zscore)) * 2


# Moderate
out.h2.women.men.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'moderate')$avg_male_DF
                                         , x = subset(d2, codeyear %in% centParties & period == 'moderate')$run_var_fem
                                         , p = 1
                                         , c = 0	
                                         , cluster = subset(d2, codeyear %in% centParties & period == 'moderate')$codeyear
                                         , all = TRUE)
out.h2.women.women.withLeg.opt <- rdrobust(y = subset(d2, codeyear %in% centParties & period == 'moderate')$avg_fem_DF
                                           , x = subset(d2, codeyear %in% centParties & period == 'moderate')$run_var_fem
                                           , p = 1
                                           , c = 0	
                                           , cluster = subset(d2, codeyear %in% centParties & period == 'moderate')$codeyear
                                           , all = TRUE)
summary(out.h2.women.men.withLeg.opt)
summary(out.h2.women.women.withLeg.opt)
# Difference 
outZ <- zDif(beta1 = out.h2.women.men.withLeg.opt$coef[3,], beta2 = out.h2.women.women.withLeg.opt$coef[3,]
             , se1 = out.h2.women.men.withLeg.opt$se[3,], se2 = out.h2.women.women.withLeg.opt$se[3,])
outZ
(1 - pnorm(outZ$zscore)) * 2

################
## Figure 4 when the mayor is a woman
################
# Get point estimates, se, and ci for men
betas_men <- c(out.h2.women.men.non.opt$coef[3,], out.h2.women.men.with.opt$coef[3,]
               , out.h2.women.men.withLeg.opt$coef[3,])
se_men <- c(out.h2.women.men.non.opt$se[3,], out.h2.women.men.with.opt$se[3,]
            , out.h2.women.men.withLeg.opt$se[3,])
lb_men <- c(out.h2.women.men.non.opt$ci[3,1], out.h2.women.men.with.opt$ci[3,1]
            , out.h2.women.men.withLeg.opt$ci[3,1])
ub_men <- c(out.h2.women.men.non.opt$ci[3,2], out.h2.women.men.with.opt$ci[3,2]
            , out.h2.women.men.withLeg.opt$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.h2.women.women.non.opt$coef[3,], out.h2.women.women.with.opt$coef[3,]
                 , out.h2.women.women.withLeg.opt$coef[3,])
se_women <- c(out.h2.women.women.non.opt$se[3,], out.h2.women.women.with.opt$se[3,]
              , out.h2.women.women.withLeg.opt$se[3,])
lb_women <- c(out.h2.women.women.non.opt$ci[3,1], out.h2.women.women.with.opt$ci[3,1]
              , out.h2.women.women.withLeg.opt$ci[3,1])
ub_women <- c(out.h2.women.women.non.opt$ci[3,2], out.h2.women.women.with.opt$ci[3,2]
              , out.h2.women.women.withLeg.opt$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.h2.women.men.non.opt$coef[3,], beta2 = out.h2.women.women.non.opt$coef[3,]
                 , se1 = out.h2.women.men.non.opt$se[3,], se2 = out.h2.women.women.non.opt$se[3,])
zOutStrong <- zDif(beta1 = out.h2.women.men.with.opt$coef[3,], beta2 = out.h2.women.women.with.opt$coef[3,]
                   , se1 = out.h2.women.men.with.opt$se[3,], se2 = out.h2.women.women.with.opt$se[3,])
zOutModerate <- zDif(beta1 = out.h2.women.men.withLeg.opt$coef[3,], beta2 = out.h2.women.women.withLeg.opt$coef[3,]
                     , se1 = out.h2.women.men.withLeg.opt$se[3,], se2 = out.h2.women.women.withLeg.opt$se[3,])

# Start the plot:
# FIGURE 4.b
pdf('./MainPaperResults/fig4_b.pdf', width = 8)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
     pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-10, 10),
     ylab = "RD Effect of Co-Partisan Mayor",
     xlab = "Electoral Rule Regime", 
     cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
                              y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
                              y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = 7.75, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
     , y = 6.5, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = 7.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
     , y = 6.5, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = 7.75, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
     , y = 6.5, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -4.5, 
       legend = c('Women', 'Men'), 
       col = c('gray46', 'black'), 
       pch = c(17, 20), 
       lty = c(1, 1),
       bty = 'n',
       cex = 1.5)
dev.off()

